#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "mathutil.h"
#include "imgutil.h"
#include "memutil.h"
#include "segment.h"
#include "global.h"

int				MM = 0;
unsigned char	*rmap9;
char			tmpfname[200];

int segment
(
	unsigned char	*rmap0,
	unsigned char	*cmap,
	int				N,
	int				ny,
	int				nx,
	int				dim,
	int				NSCALE
)
{
	/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
	float	MINRSIZE[5];
	int		scale[5], offset[5], step[5], MAXSCALE, MINSCALE, *count, TR, i, j, k, l;
	int		oldTR, extraTR, *reg, *reg2, *convert, datasize, autoscale;
	char	fname[200];
	short	*rmap;
	/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/

	
	step[0] = 1;	// not used
	step[1] = 1;
	step[2] = 2;
	step[3] = 4;
	step[4] = 8;

	scale[0] = 32;
	scale[1] = 64;	
	scale[2] = 128;
	scale[3] = 256;
	scale[4] = 512;
	
	offset[0] = 2;	
	offset[1] = 4;
	offset[2] = 8;	
	offset[3] = 16;	
	offset[4] = 32;
		
	for(i = 4; i >= 1; i--)
	{
		if(ny * nx >= sqr(scale[i]))
		{
			MAXSCALE = i;
			break;
		}
	}

	if(i == 0)
	{
		return 0;
	}

	for(i = 0; i <= MAXSCALE; i++)
	{
		MINRSIZE[i] = 2.0 * sqr(offset[i]);
	}

	if(NSCALE > 0)
	{
		autoscale = 0;
		MINSCALE = MAXSCALE - NSCALE + 1;
		if(MINSCALE <= 0)
		{
			MINSCALE = 1;
			NSCALE = MAXSCALE - MINSCALE + 1;
		}
	}
	else if(NSCALE == -1)
	{
		autoscale = 1;
		if(MAXSCALE >= 2)
		{
			MINSCALE = 2;
			NSCALE = MAXSCALE - MINSCALE + 1;
			if(NSCALE > 2)
			{
				NSCALE = 2;
				MINSCALE = MAXSCALE - NSCALE + 1;
			}
		}
		else
		{
			MINSCALE = 1;
			NSCALE = 1;
		}
	}

	datasize = ny * nx;
	rmap = (short *) calloc(datasize, sizeof(short));
	for(l = 0; l < datasize; l++)
	{
		rmap0[l] = 1;
	}

	oldTR = 1;

	for(i = MAXSCALE; i >= MINSCALE; i--)
	{
		for(l = 0; l < datasize; l++)
		{
			rmap[l] = 0;
		}

		TR = segment1(cmap, N, ny, nx, offset, step, rmap, rmap0, oldTR, i, MINRSIZE[i]);

		reg = (int *) calloc(oldTR + 1, sizeof(int));
		reg2 = (int *) calloc(TR + 1, sizeof(int));
		for(l = 0; l < datasize; l++)
		{
			j = rmap0[l];
			k = rmap[l];
			if(reg[j] == 0)
			{
				reg[j] = k;
			}
			else if(reg[j] == -1)
			{
				reg2[k] = -1;
			}
			else if(reg[j] != k)
			{
				reg2[k] = -1;
				reg2[reg[j]] = -1;
				reg[j] = -1;
			}

			rmap0[l] = k;
		}

		free(reg);

		count = (int *) calloc(TR + 1, sizeof(int));
		for(l = 0; l < datasize; l++)
		{
			count[rmap0[l]]++;
		}

		do
		{
			if(i <= MAXSCALE - 2 || i == 1)
			{
				break;
			}

			oldTR = TR;
			convert = (int *) calloc(oldTR + 1, sizeof(int));
			for(j = 1; j <= oldTR; j++)
			{
				convert[j] = j;
			}

			for(j = 1; j <= oldTR; j++)
			{
				if(count[j] > sqr(scale[i]) / 8 && reg2[j] == -1)
				{
					convert[j] = 0;
					TR--;
					for(k = j + 1; k <= oldTR; k++)
					{
						convert[k]--;
					}
				}
			}

			if(TR < oldTR)
			{
				for(l = 0; l < datasize; l++)
				{
					if(convert[rmap0[l]] == 0)
					{
						rmap[l] = 0;
					}
					else
					{
						rmap[l] = -1;
					}
				}

				extraTR = segment1(cmap, N, ny, nx, offset, step, rmap, rmap0, oldTR, i, MINRSIZE[i]);
				if(extraTR > oldTR - TR)
				{
					reg2 = (int *) realloc(reg2, (TR + extraTR + 1) * sizeof(int));
					reg = (int *) calloc(oldTR + 1, sizeof(int));
					for(j = 1; j <= oldTR; j++)
					{
						reg2[convert[j]] = reg2[j];
					}

					for(j = TR + 1; j <= TR + extraTR; j++)
					{
						reg2[j] = 0;
					}

					for(l = 0; l < datasize; l++)
					{
						j = rmap0[l];
						if(rmap[l] > 0)
						{
							k = TR + rmap[l];
							rmap0[l] = k;
							if(reg[j] == 0)
							{
								reg[j] = k;
							}
							else if(reg[j] == -1)
							{
								reg2[k] = -1;
							}
							else if(reg[j] != k)
							{
								reg2[k] = -1;
								reg2[reg[j]] = -1;
								reg[j] = -1;
							}
						}
						else if(rmap[l] == 0)
						{
							rmap0[l] = 0;
						}
						else
						{
							rmap0[l] = convert[j];
						}
					}

					free(reg);
				}

				TR += extraTR;
			}

			free(convert);

			if(TR > oldTR)
			{
				count = (int *) realloc(count, (TR + 1) * sizeof(int));
				for(j = 1; j <= TR; j++)
				{
					count[j] = 0;
				}

				for(l = 0; l < datasize; l++)
				{
					count[rmap0[l]]++;
				}
			}

		} while(oldTR != TR);
		free(count);
		free(reg2);

		oldTR = TR;
	}

	if(autoscale == 1 && MINSCALE > 1 && NSCALE < 3)
	{
		i = MINSCALE - 1;
		TR = segment2(rmap0, rmap, i, cmap, N, ny, nx, oldTR, MINRSIZE, offset, step);
	}

	free(rmap);
	return TR;
}

int segment2
(
	unsigned char	*rmap0,
	short			*rmap,
	int				i,
	unsigned char	*cmap,
	int				N,
	int				ny,
	int				nx,
	int				oldTR,
	float			*MINRSIZE,
	int				*offset,
	int				*step
)
{
	/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
	float	oldJ[256], **mapmatrix, *change, **P;
	int		TR, j, k, l, *convert, imgsize, *count, extraTR;
	/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
	
	imgsize = ny * nx;
	for(k = 0; k < 256; k++)
	{
		oldJ[k] = 0;
	}
	
	mapmatrix = (float **) fmatrix(oldTR + 1, N + 1);
	gettotalJS(cmap, N, ny, nx, rmap0, oldTR, oldJ, mapmatrix, oldTR);

	count = (int *) calloc(oldTR + 1, sizeof(int));
	P = (float **) fmatrix(oldTR + 1, N);
	for(l = 0; l < imgsize; l++)
	{
		P[rmap0[l]][cmap[l]] += 1.0;
		count[rmap0[l]]++;
	}

	for(j = 1; j <= oldTR; j++)
	{
		for(k = 0; k < N; k++)
		{
			mapmatrix[j][k] /= P[j][k];
		}

		mapmatrix[j][N] /= count[j];
	}

	change = (float *) calloc(oldTR + 1, sizeof(float));
	convert = (int *) calloc(oldTR + 1, sizeof(int));
	for(j = 1; j <= oldTR; j++)
	{
		for(k = 0; k < N; k++)
		{
			if(mapmatrix[j][k] < mapmatrix[j][N] / 4 && P[j][k] > 2 * MINRSIZE[i])
			{
				change[j] = 1.0;
			}
		}
	}

	TR = oldTR;
	for(j = 1; j <= oldTR; j++)
	{
		convert[j] = j;
	}

	for(j = 1; j <= oldTR; j++)
	{
		if(change[j] >= 1.0 || oldJ[j] > 1.0)
		{
			convert[j] = 0;
			for(k = j + 1; k <= oldTR; k++)
			{
				convert[k]--;
			}

			TR--;
		}
	}

	for(l = 0; l < imgsize; l++)
	{
		if(convert[rmap0[l]] > 0)
		{
			rmap[l] = -1;
		}
		else
		{
			rmap[l] = 0;
		}
	}

	extraTR = segment1(cmap, N, ny, nx, offset, step, rmap, rmap0, oldTR, i, MINRSIZE[i]);
	for(l = 0; l < imgsize; l++)
	{
		if(rmap[l] > 0)
		{
			rmap0[l] = rmap[l] + TR;
		}
		else
		{
			rmap0[l] = convert[rmap0[l]];
		}
	}

	TR += extraTR;

	free(count);
	free(convert);
	free(change);
	free_fmatrix(P, oldTR + 1);
	free_fmatrix(mapmatrix, oldTR + 1);
	return TR;
}

int segment1
(
	unsigned char	*cmap,
	int				N,
	int				ny,
	int				nx,
	int				*offset,
	int				*step,
	short			*rmap,
	unsigned char	*rmap0,
	int				oldTR,
	int				i,
	float			MINRSIZE
)
{
	/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
	float	*J, *JT, *threshJ1, *threshJ2;
	int		j, TR, *appear, *done, alldone, datasize, imgsize, l;
	/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/

	imgsize = ny * nx;
	datasize = ny * nx;
	J = (float *) calloc(datasize, sizeof(float));

	getJ
	(
		cmap,
		N,
		ny,
		nx,
		J,
		offset[i],
		step[i],
		rmap,
		rmap0,
		oldTR
	);
	
	JT = (float *) calloc(0, sizeof(float));

	done = (int *) calloc(oldTR + 1, sizeof(int));
	threshJ1 = (float *) calloc(oldTR + 1, sizeof(float));
	threshJ2 = (float *) calloc(oldTR + 1, sizeof(float));
	
	getthreshJ
	(
		imgsize,
		J,
		rmap,
		rmap0,
		threshJ1,
		threshJ2,
		oldTR,
		0,
		done
	);
	
	appear = (int *) calloc(oldTR + 1, sizeof(int));
	TR = getrmap3(rmap, J, ny, nx, threshJ1, 0, MINRSIZE, rmap0, 0, threshJ2, oldTR, appear);
	free(appear);

	free(threshJ1);
	free(threshJ2);
	free(JT);

	removehole(rmap, ny, nx, rmap0);
	alldone = getrmap2(rmap, J, ny, nx, TR, oldTR, rmap0, done);
	if(alldone < oldTR)
	{
		removehole(rmap, ny, nx, rmap0);
		for(j = i - 1; j >= MAX(i - 2, 1); j--)
		{
			getJ
			(
				cmap,
				N,
				ny,
				nx,
				J,
				offset[j],
				step[j],
				rmap,
				rmap0,
				oldTR
			);

			alldone = getrmap2(rmap, J, ny, nx, TR, oldTR, rmap0, done);
			if(alldone == oldTR)
			{
				break;
			}

			removehole(rmap, ny, nx, rmap0);
		}

		flood(rmap, J, ny, nx, rmap0, oldTR, done);
	}

	free(done);
	free(J);

	return TR;
}


int merge1(unsigned char *rmap, unsigned char *cmap, int N, int ny, int nx, int TR, float threshcolor)
{
	/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
	int				iy, ix, ir, jr, i, mini, minj, newtr, *npt, npttotal;
	float			*currentJ, *mergeJ, mindist, **distnpt, **distcolor, **distJ, **unused, **P;
	int				loc, loc1, l, datasize, imgsize, *convert, threshtr;
	unsigned char	*rmap1, *rmap2;
	float			oldoverallJ, overallJ;
	/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/

	if(threshcolor >= 0)
	{
		TR = merge(rmap, cmap, N, ny, nx, TR, threshcolor, 0);
		return TR;
	}

	threshcolor = 0.5;
	threshtr = TR;

	imgsize = ny * nx;
	datasize = imgsize;
	distnpt = (float **) fmatrix(TR + 1, TR + 1);
	distJ = (float **) fmatrix(TR + 1, TR + 1);
	distcolor = (float **) fmatrix(TR + 1, TR + 1);
	loc = 0;
	
	for(iy = 0; iy < ny; iy++)
	{
		loc++;
		for(ix = 1; ix < nx; ix++)
		{
			loc1 = loc - 1;
			if(rmap[loc1] != rmap[loc])
			{
				ir = rmap[loc1];
				jr = rmap[loc];
				distnpt[ir][jr] += 1;
				distnpt[jr][ir] = distnpt[ir][jr];
			}

			loc++;
		}
	}

	loc = 0;
	
	loc += nx;
	for(iy = 1; iy < ny; iy++)
	{
		for(ix = 0; ix < nx; ix++)
		{
			loc1 = loc - nx;
			if(rmap[loc1] != rmap[loc])
			{
				ir = rmap[loc1];
				jr = rmap[loc];
				distnpt[ir][jr] += 1;
				distnpt[jr][ir] = distnpt[ir][jr];
			}

			loc++;
		}
	}

	rmap1 = (unsigned char *) calloc(datasize, sizeof(unsigned char));
	currentJ = (float *) calloc(TR + 1, sizeof(float));
	mergeJ = (float *) calloc(2, sizeof(float));
	overallJ = gettotalJS(cmap, N, ny, nx, rmap, TR, currentJ, unused, 0);
	oldoverallJ = overallJ;

	P = (float **) fmatrix(TR + 1, N);
	npt = (int *) calloc(TR + 1, sizeof(int));
	for(l = 0; l < datasize; l++)
	{
		P[rmap[l]][cmap[l]] += 1;
		npt[rmap[l]]++;
	}

	for(ir = 1; ir <= TR; ir++)
	{
		for(i = 0; i < N; i++)
		{
			P[ir][i] /= npt[ir];
		}
	}

	mindist = 10000;
	for(ir = 1; ir < TR; ir++)
	{
		for(jr = ir + 1; jr <= TR; jr++)
		{
			if(distnpt[ir][jr] != 0.0)
			{
				for(l = 0; l < datasize; l++)
				{
					if(rmap[l] == ir || rmap[l] == jr)
					{
						rmap1[l] = 1;
					}
					else
					{
						rmap1[l] = 0;
					}
				}

				gettotalJS(cmap, N, ny, nx, rmap1, 1, mergeJ, unused, 0);
				distJ[ir][jr] = mergeJ[1] - (npt[ir] * currentJ[ir] + npt[jr] * currentJ[jr]) / (npt[ir] + npt[jr]);
				distJ[jr][ir] = distJ[ir][jr];

				distcolor[jr][ir] = distcolor[ir][jr] = distance(P[ir], P[jr], N);
				if(distcolor[ir][jr] < mindist)
				{
					mindist = distcolor[ir][jr];
					mini = ir;
					minj = jr;
				}
			}
		}
	}

	convert = (int *) calloc(TR + 1, sizeof(int));
	for(i = 1; i <= TR; i++)
	{
		convert[i] = i;
	}

	newtr = TR;
	while(mindist < threshcolor)
	{
		overallJ = overallJ - (npt[mini] * currentJ[mini] + npt[minj] * currentJ[minj]) / datasize;
		npttotal = npt[mini] + npt[minj];
		for(i = 0; i < N; i++)
		{
			P[mini][i] = (P[mini][i] * npt[mini] + P[minj][i] * npt[minj]) / npttotal;
		}

		currentJ[mini] = (npt[mini] * currentJ[mini] + npt[minj] * currentJ[minj]) / npttotal + distJ[mini][minj];
		npt[mini] = npttotal;
		overallJ = overallJ + npt[mini] * currentJ[mini] / datasize;

		if(overallJ <= oldoverallJ)
		{
			oldoverallJ = overallJ;
			threshtr = newtr - 1;
		}

		for(i = 1; i <= TR; i++)
		{
			if(convert[i] == minj)
			{
				convert[i] = mini;
			}
		}

		for(ir = 1; ir <= newtr; ir++)
		{
			if(ir != mini && ir != minj)
			{
				if(distnpt[mini][ir] != 0.0 || distnpt[minj][ir] != 0.0)
				{
					for(l = 0; l < datasize; l++)
					{
						if(convert[rmap[l]] == mini || convert[rmap[l]] == ir)
						{
							rmap1[l] = 1;
						}
						else
						{
							rmap1[l] = 0;
						}
					}

					gettotalJS(cmap, N, ny, nx, rmap1, 1, mergeJ, unused, 0);
					distJ[mini][ir] = mergeJ[1] - (npt[mini] * currentJ[mini] + npt[ir] * currentJ[ir]) / (npt[mini] + npt[ir]);
					distJ[ir][mini] = distJ[mini][ir];

					distcolor[mini][ir] = distance(P[mini], P[ir], N);
					distcolor[ir][mini] = distcolor[mini][ir];
					distnpt[mini][ir] = distnpt[mini][ir] + distnpt[minj][ir];
					distnpt[ir][mini] = distnpt[mini][ir];
				}
			}
		}

		for(i = 1; i <= TR; i++)
		{
			if(convert[i] > minj)
			{
				convert[i]--;
			}
		}

		for(ir = minj + 1; ir <= newtr; ir++)
		{
			npt[ir - 1] = npt[ir];
			currentJ[ir - 1] = currentJ[ir];
			for(i = 0; i < N; i++)
			{
				P[ir - 1][i] = P[ir][i];
			}
		}

		for(ir = 1; ir <= newtr; ir++)
		{
			for(jr = minj + 1; jr <= newtr; jr++)
			{
				distnpt[ir][jr - 1] = distnpt[ir][jr];
				distcolor[ir][jr - 1] = distcolor[ir][jr];
				distJ[ir][jr - 1] = distJ[ir][jr];
			}
		}

		for(ir = minj + 1; ir <= newtr; ir++)
		{
			for(jr = 1; jr <= newtr - 1; jr++)
			{
				distnpt[ir - 1][jr] = distnpt[ir][jr];
				distcolor[ir - 1][jr] = distcolor[ir][jr];
				distJ[ir - 1][jr] = distJ[ir][jr];
			}
		}

		newtr--;
		mindist = 10000;

		for(ir = 1; ir < newtr; ir++)
		{
			for(jr = ir + 1; jr <= newtr; jr++)
			{
				if(distnpt[ir][jr] != 0)
				{
					if(distcolor[ir][jr] < mindist)
					{
						mindist = distcolor[ir][jr];
						mini = ir;
						minj = jr;
					}
				}
			}
		}
	}

	/*
	 * for (l=0;
	 * l<datasize;
	 * l++) rmap[l]=convert[rmap[l]];
	 */
	free(rmap1);
	free(currentJ);
	free(mergeJ);
	free(npt);
	free_fmatrix(P, TR + 1);
	free(convert);
	free_fmatrix(distnpt, TR + 1);
	free_fmatrix(distJ, TR + 1);
	free_fmatrix(distcolor, TR + 1);

	TR = merge(rmap, cmap, N, ny, nx, TR, threshcolor, threshtr);
	return TR;
}

int merge
(
	unsigned char	*rmap,
	unsigned char	*cmap,
	int				N,
	int				ny,
	int				nx,
	int				TR,
	float			threshcolor,
	int				threshtr
)
{
	/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
	int		iy, ix, ir, *npt, npttotal, jr, i, mini, minj, newtr;
	float	**P, mindist, **distnpt, **distcolor, nptT = 0.2;
	int		loc, loc1, l, datasize, imgsize, *convert;;
	/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/

	imgsize = ny * nx;
	datasize = imgsize;
	distnpt = (float **) fmatrix(TR + 1, TR + 1);
	distcolor = (float **) fmatrix(TR + 1, TR + 1);
	loc = 0;
	
	for(iy = 0; iy < ny; iy++)
	{
		loc++;
		for(ix = 1; ix < nx; ix++)
		{
			loc1 = loc - 1;
			if(rmap[loc1] != rmap[loc])
			{
				ir = rmap[loc1];
				jr = rmap[loc];
				distnpt[ir][jr] += 1;
				distnpt[jr][ir] = distnpt[ir][jr];
			}

			loc++;
		}
	}

	loc = 0;
	
	loc += nx;
	for(iy = 1; iy < ny; iy++)
	{
		for(ix = 0; ix < nx; ix++)
		{
			loc1 = loc - nx;
			if(rmap[loc1] != rmap[loc])
			{
				ir = rmap[loc1];
				jr = rmap[loc];
				distnpt[ir][jr] += 1;
				distnpt[jr][ir] = distnpt[ir][jr];
			}

			loc++;
		}
	}

	P = (float **) fmatrix(TR + 1, N);
	npt = (int *) calloc(TR + 1, sizeof(int));
	for(l = 0; l < datasize; l++)
	{
		ir = rmap[l];
		P[ir][cmap[l]] += 1;
		npt[ir]++;
	}

	for(ir = 1; ir <= TR; ir++)
	{
		for(i = 0; i < N; i++)
		{
			P[ir][i] /= npt[ir];
		}
	}

	mindist = 10000;
	for(ir = 1; ir < TR; ir++)
	{
		for(jr = ir + 1; jr <= TR; jr++)
		{
			if(distnpt[ir][jr] != 0.0)
			{
				distcolor[jr][ir] = distcolor[ir][jr] = distance(P[ir], P[jr], N);
				if(distcolor[ir][jr] < mindist)
				{
					mindist = distcolor[ir][jr];
					mini = ir;
					minj = jr;
				}
			}
		}
	}

	convert = (int *) calloc(TR + 1, sizeof(int));
	for(i = 1; i <= TR; i++)
	{
		convert[i] = i;
	}

	newtr = TR;
	while(mindist < threshcolor && newtr > threshtr)
	{
		for(i = 1; i <= TR; i++)
		{
			if(convert[i] == minj)
			{
				convert[i] = mini;
			}
			else if(convert[i] > minj)
			{
				convert[i]--;
			}
		}

		npttotal = npt[mini] + npt[minj];
		for(i = 0; i < N; i++)
		{
			P[mini][i] = (P[mini][i] * npt[mini] + P[minj][i] * npt[minj]) / npttotal;
		}

		npt[mini] = npttotal;
		for(ir = 1; ir <= newtr; ir++)
		{
			if(ir != mini && ir != minj)
			{
				if(distnpt[mini][ir] != 0.0 || distnpt[minj][ir] != 0.0)
				{
					distcolor[mini][ir] = distance(P[mini], P[ir], N);
					distcolor[ir][mini] = distcolor[mini][ir];
					distnpt[mini][ir] = distnpt[mini][ir] + distnpt[minj][ir];
					distnpt[ir][mini] = distnpt[mini][ir];
				}
			}
		}

		for(ir = 1; ir <= newtr; ir++)
		{
			for(jr = minj + 1; jr <= newtr; jr++)
			{
				distnpt[ir][jr - 1] = distnpt[ir][jr];
				distcolor[ir][jr - 1] = distcolor[ir][jr];
			}
		}

		for(ir = minj + 1; ir <= newtr; ir++)
		{
			for(jr = 1; jr <= newtr - 1; jr++)
			{
				distnpt[ir - 1][jr] = distnpt[ir][jr];
				distcolor[ir - 1][jr] = distcolor[ir][jr];
			}

			npt[ir - 1] = npt[ir];
			for(i = 0; i < N; i++)
			{
				P[ir - 1][i] = P[ir][i];
			}
		}

		newtr--;
		mindist = 10000;
		for(ir = 1; ir < newtr; ir++)
		{
			for(jr = ir + 1; jr <= newtr; jr++)
			{
				if(distnpt[ir][jr] != 0)
				{
					if(distcolor[ir][jr] < mindist)
					{
						mindist = distcolor[ir][jr];
						mini = ir;
						minj = jr;
					}
				}
			}
		}
	}

	for(l = 0; l < datasize; l++)
	{
		rmap[l] = convert[rmap[l]];
	}

	free(convert);
	free(npt);
	free_fmatrix(P, TR + 1);
	free_fmatrix(distnpt, TR + 1);
	free_fmatrix(distcolor, TR + 1);
	return newtr;
}
